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We investigate by direct numerical solution of appropriate renormalization flow equations, the 
validity of a recent dislocation unbinding theory for laser induced freezing/melting in two dimensions. 
The bare elastic moduli and dislocation fugacities which are inputs to the flow equations are obtained 
for three different 2-d systems (hard disk, inverse 12 th power and the Derjaguin-Landau-Verwey- 
Overbeek potentials) from a restricted Monte Carlo simulation sampling only configurations without 
dislocations. We conclude that (a) the flow equations need to be correct at least up to third order 
in defect fugacity to reproduce meaningful results, (b) there is excellent quantitative agreement 
between our results and earlier conventional Monte Carlo simulations for the hard disk system and 
(c) while the qualitative form of the phase diagram is reproduced for systems with soft potentials 
there is some quantitative discrepancy which we explain. 
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I. INTRODUCTION 

Examples of phase transitions mediated by the unbind- 
ing of defect pairs abound in two dimensions. The quasi- 
longranged- order to disorder transition in the XY and 
planar rotor models^, the melting transition of a two 
dimensional solid 3 -, the superconductor to normal phase 
transition in two dimensional Joscphson junction arrays*, 
the commensurate- incommensurate transition of the 
striped phase of smectic liquid crystals on anisotropic 
substrates 5 -, and the more recent discovery of a defect me- 
diated re-entrant freezing transition in two dimensional 
colloids in an external periodic potential^ are all under- 
stood within such defect unbinding theories. While the 
very first defect mediated transition theory for the phase 
transition in the XY-model by Kosterlitz and Thouless 
(KT)i enjoyed almost immediate acceptance and was 
verified in simulations 2 as well as experiments^, de- 
fect mediated theories of two dimensional melting took a 
long time to gain general acceptance in the community 10 . 
There were several valid reasons for this reticence how- 
ever. 

Firstly, as was recognized even in the earliest papersii^ 
on this subject, the dislocation unbinding transition, 
which represents an instability of the solid phase, may al- 
ways be pre-empted by a first ordeni&ia transition from 
a metastable solid to a stable liquid. Whether such a 
first order melting transition actually occurs or not de- 
pends on the temperature of instability Tkt] so that if 
the transition temperature T c < Tkt the unbinding of 
dislocations does not occur. Clearly, neither this condi- 
tion nor its converse can hold for all 2d systems in gen- 
eral since Tkt is a non-universal number which depends 
on the "distance" in coupling parameter space between 
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the bare and the fixed point Hamiltonian and hence on 
the details of the interaction. Secondly, the renormaliza- 
tion group flow equations derived in all defect mediated 
theories to date are perturbative expansions in the de- 
fect density (fugacity) in the ordered phase. How fast 
does this perturbation series converge? Again, the an- 
swer lies in the position of the bare Hamiltonian in the 
coupling parameter space. For the planar rotor modeU^, 
past calculations show that next to leading order terms 
in the flow equations are essential to reproduce the value 
of the transition temperature obtained in simulations 2 . 
Thirdly, defect mediated transitions predict an essen- 
tial singularity^ of the correlation length at the transi- 
tion temperature. This means that effects of finite size 15 
would be substantial and may thoroughly mask the true 
thermodynamic result. A rapid increase of the correla- 
tion length also implies that the relaxation time diverges 
as the transition temperature is approached - critical 
slowing down. For the two dimensional solid, this last 
effect is particularly crucial since, even far from the tran- 
sition, the motion of defects is mainly thermally assisted 
and diffusional and therefore slow. The equilibration of 
defect configurations— is therefore often an issue even in 
solids of macroscopic dimensions. 

On the other hand, over the last few years it has been 
possible to test quantitatively some of the non-universal 
predictions of defect mediated theories of phase tran- 
sitions using simulations of restricted systems 2 * 3 *^. A 
simulation of a system without defects is used to ob- 
tain the values for the bare coupling constants which are 
then taken as inputs to the renormalization group equa- 
tions for the appropriate defect unbinding theory to ob- 
tain quantities like the transition temperature. Needless 
to say, the simulated system does not undergo a phase 
transition and therefore problems typically related to di- 
verging correlation lengths and times do not occur. Nu- 
merical agreement of the result of this calculation with 
that of unrestricted simulations or experiments is proof 
of the validity of the RG flow equations 1 : 3 : 11 ! 12 . This 
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idea has been repeatedly applied in the past to ana- 
lyze defect mediated phase transitions in the planar rotor 
model^, two dimensional melting of hard disks?» and the 
re-entrant freezing of hard disks in an external periodic 
potentialiLi£. The last system is particularly interesting 
in view of its close relation with experiments on laser in- 
duced re-entrant freezing transition in charge stabilized 
colloids^ and this constitutes the subject of the present 
paper as well. 

In this paper we show in detail how restricted simulations 
of systems of particles interacting among themselves via a 
variety of interactions and with a commensurate external 
periodic potential can be used to obtain phase diagrams 
showing the re-entrant freezing transition. The results 
obtained are compared to earlier unrestricted simulations 
for the same systems. Briefly our results are as follows. 
Firstly, we observe that, as in an earlier study of the pla- 
nar rotor model^, next to leading order corrections to the 
renormalization flow equations are essential to reproduce 
even the gross features of the phase diagram. Specifi- 
cally, the re-entrant portion of the phase diagram can be 
reproduced only if such correction terms are taken into 
account. Secondly while we find almost complete agree- 
ment with earlier results for the hard disk system which 
has been studied most extensively, our phase diagram for 
the other forms of interaction is shifted with respect to 
the results available in the literature. This may mean ei- 
ther of two things — inadequacy of the RG theory used 
by us or finite size effects in the earlier results. Lastly, as 
a by product of our calculations, we have obtained the 
core energy for defects (dislocations) in these systems and 
studied its dependence on thermodynamic and potential 
parameters. 

The problem of re-entrant freezing transition of a sys- 
tem of interacting colloidal particles in a periodic poten- 
tial has an interesting history involving experiments^ 7 -, 
simulation 1 V"i 21 i 22 i 23 i 24 i 25 i 26 and theory^. In last 
couple of decades soft systems like colloids have been 
studied extensively 2 ^ both for their own sake and as 
typical toy models to study various important con- 
densed matter questions like structural and phase transi- 
tions through experiments that allow real space imaging. 
Charged colloids confined within two glass plates form 
a model 2-d system as the electrostatic force from the 
plates almost completely suppresses the fluctuations of 
colloids perpendicular to the plates, practically confin- 
ing them to a 2-d plane. In their pioneering experiment 
Chowdhury 6 et. al. imposed a simple static background 
potential which is periodic in one direction and constant 
in the other (except for an overall Gaussian profile of 
intensity- variation) by interfering two laser beams. This 
potential immediately induces a density modulation in 
the colloidal system. The potential minima are spaced 
to overlap with the close packed lines of the ideal lattice 
of the colloidal system at a given density. With increase 
in potential strength such a colloidal liquid has been ob- 
served to solidify. This is known as laser induced freezing 
(LIF). In a recent experiment! it has been shown that 



with further increase in potential strength, surprisingly, 
the solid phase re-melts to a modulated liquid. This phe- 
nomenon is known as re-entrant laser induced freezing 
(RLIF). Qualitatively, starting from a liquid phase, the 
external periodic potential immediately induces a density 
modulation, reducing fluctuations which eventually leads 
to solidification. Further increase in the amplitude of the 
potential reduces the system to a collection of decoupled 
1-d strips. The dimensional reduction now increases fluc- 
tuations remelting the system. 

The early mean field theories, namely, Landau theoryS 
and density functional theory 2 ^ predicted a change from 
a first order to continuous transition with increase in 
potential strength and failed to describe the re-entrant 
behavior, a conclusion seemingly confirmed by early 
experiments^ and some early simulations^. Overall, the 
results from early simulations remained inconclusive how- 
ever, while one of them 19 claimed to have found a tri- 
critical point at intermediate laser intensities and re- 
entrance, later studies refuted these result o 20 ' 21 i 22 . All 
of these studies used the change in order parameter and 
the maximum in the specific heat to identify the phase 
transition points. While the later studiesS&SiiS^ found 
RLIF for hard disks they reported laser induced freez- 
ing and absence of any re-entrant melting for the DLVO 
potential 22 in direct contradiction to experiments 7 . 

Following the defect mediated disordering approach 
of Kosterlitz and Thouless 1 (KT), Frey, Nelson and 
Radzihovsky^S (FNR) proposed a detailed theory for the 
re-entrant transition based on the unbinding of disloca- 
tions with Burger's vector parallel to the line of potential 
minima. This theory predicted RLIF and no tricritical 
point. The results of this work were in qualitative agree- 
ment with experiments- and provided a framework for 
understanding RLIF in general. More accurate simula- 
tion studies on systems of hard disks^ 3 -, soft disks^Si, 
DLVO 2 - etc. confirmed the re-entrant freezing-melting 
transition in agreement with experiments* 7 - and FNR 
theory 2 ^. In these studies the phase transition point was 
found from the crossing of Binder-cumulants^S^i of order 
parameters corresponding to translational and bond- ori- 
entational order, calculated for various sub- system sizes. 
A systematic finite size scaling analysis 2 ^ of simulation 
results for the 2-d hard disk system in a 1-d modulat- 
ing potential showed, in fact, several universal features 
consistent with the predictions of FNR theory. It was 
shown in these studies that the resultant phase diagram 
remains system size dependent and the cross- over to 
the zero field KTHNY meltingiiiiS plays a crucial role 
in understanding the results for small values of the ex- 
ternal potential. While the data collapse and critical 
exponents were consistent with KT theory for stronger 
potentials, for weaker potentials they match better with 
critical scaling 23 . A common problem with all the simula- 
tion studies might be equilibration with respect to dislo- 
cation movements along climb (or even glide) directions. 
Also, non universal predictions, namely the phase dia- 
gram are difficult to compare because the FNR approach 
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FIG. 1: This cartoon shows a typical 2-d system under consid- 
eration, d is the length scale over which repulsive two body 
potentials are operative. The dashed lines indicate minima 
of external modulating potential (3V(y) = —/3Vo cos(2ny / do) ■ 
a x — ao is the lattice parameter fixed by the density p and 
a y indicate the average separation between two layers along 
y-direction perpendicular to a set of close-packed planes. For 
a perfect triangular lattice a y — v3<Jo/2. The modulating 
potential is commensurate with the lattice such that do — a y . 



(like KT theory) is expressed in terms of the appropriate 
clastic moduli which are notoriously time-consuming to 
compute near a continuous phase transition. Diverging 
correlation lengths and times near the phase transition 
point further complicate an accurate evaluation of the 
non universal predictions of the theory. 

We calculate the phase diagrams of three different 
2-d systems with a 1-d modulating potential (see 
Fig. ^| following the technique of restricted Monte Carlo 
simulations^iii, to be discussed below. For the laser in- 
duced transition we use this method to generate whole 
phase diagrams. We reject Monte Carlo moves which 
tend to distort an unit cell in a way which changes the 
local connectivity^. The percentage of moves thus re- 
jected is a measure of the dislocation fugacitj^. This, 
together with the elastic constants of the dislocation free 
lattice obtained separately, are our inputs (bare values) 
to the renormalization flow equations 28 to compute the 
melting points and hence the phase diagram. Our results 
(Fig. 110112111(1 clearly show a modulated liquid (ML) -> 
locked floating solid (LFS) — > ML re-entrant transition 
with increase in the amplitude (Vo) of the potential. In 
general, we find, the predictions of FNR theory to be 
valid. 

In section ITTlwc first briefly discuss the FNR theory and 
then go on to show in detail the restricted simulation 
scheme used by us to obtain the various quantities re- 
quired to calculate the phase diagram. In section ITTT1 we 
give the simulation results. We describe, in detail, the 
various quantities leading to the phase diagram for one 
of the systems, viz. the hard disks^^. Then we present 
the phase diagrams for the other two systems we study. 
We compare our results with earlier simulations. Lastly, 
in section IIVI we summarize our main results and con- 
clude. 



II. METHOD 

A cartoon corresponding to the systems considered for 
our study is given in Fig. The elastic free energy 
of the solid is given in terms of the spatial derivatives 
of the displacement field u(r) = r — r$ with being 
the lattice vectors of the undistorted reference triangular 
lattice. For a solid in presence of a modulating poten- 
tial (3V{y) (Fig. nj the displacement mode u y becomes 
massive, leaving massless u x modes. After integrating 
out the u v modes the free energy of the LFS may be ex- 
pressed in terms of gradients of u x and elastic moduli, 
namely, the Young's modulus K(f3Vo,p) and shear mod- 
ulus fi(j3Vo,p), 
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Similar arguments^ show that among the three sets of 
low energy dislocations available in the 2-d triangular lat- 
tice, only those (type I) with Burger's vector parallel to 
the line of potential minima survive at large /?Vb- Dislo- 
cations with Burger's vector pointing along the other two 
possible close-packed directions (type II) in the 2-d trian- 
gular lattice have larger energies because the surround- 
ing atoms are forced to ride the crests of the periodic 
potential 2 ^. Within this set of assumptions, the system 
therefore shares the same symmetries as the XY model. 
Indeed, a simple rescaling of x — * y/JIx and y — ► y/~Ky 
leads this free energy to the free energy of the XY-model 
with spin- wave stiffness K xy — ^/ K fia^ / 'Att 2 and spin an- 
gle 6 = 2wu x /a Q : 



dxdy 
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This immediately leads to the identification of a vortex 
in XY model (§ d9 — 2tt) with a dislocation of Burger's 

vector b = iao du x = ao, i = unit vector along x- direc- 
tion) parallel to the potential minima i.e. the dislocation 
of type I. The corresponding theory for phase transitions 
can then be recast as a KT theorji and is described in 
the framework of a two parameter renormalization flow 
for the spin- wave stiffness K xy (l) and the fugacity of type 
I dislocations y'(l), where I is a measure of length scale 
as I = ln(r/ao), r being the size of the system. The 
flow equations are expressed in terms of x' = (TrK xy — 2) 
and y' — Ait exp(—/3E c ) where E c is the core energy of 
type I dislocations which is obtained from the disloca- 
tion probability&iS. Keeping upto next to leading order 
terms in y' the renormalization group flow equations 2 ^^ 
are, 
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Flows in / generated by the above equations starting from 
initial or "bare" values of x' and y' fall in two categories. 
If, as I — > oo, y' diverges, the thermodynamic phase is 
disordered (i.e. ML), while on the other hand if y' van- 
ishes, it is an ordered phase (a LFS} 2 -. The two kinds 
of flows are demarcated by the separatrix which marks 
the phase transition point. For the linearized equations, 
that keeps upto only the leading order terms in y', the 
separatrix is simply the straight line y' = x' , whereas for 
the full non-linear equations one needs to calculate this 
numericallyVi 34 . 

The bare numbers for x' and y' are relatively insensitive 
to system size since our Monte Carlo simulation does not 
involve a diverging correlation length associated with a 
phase transition. This is achieved as follows^ 3 -. We mon- 
itor individual random moves of the particles in a system 
and look for distortions of the neighboring unit cells. If 
in any of these unit cells the length of a next nearest 
neighbor bond becomes smaller than the nearest neigh- 
bor bond, the move is rejected. All such moves generate 
disclination quartets and are shown in the Fig. [21 Notice 
that each of these moves break a nearest neighbour bond 
to build a new next nearest neighbour bond, in the pro- 
cess generating two 7-5 disclination pairs. These are the 
moves rejected in the restricted simulation scheme we 
follow. The probabilities of such bond breaking moves 
are however computed by keeping track of the number 
of such rejected moves. One has to keep track of all 
the three possible distortions of the unit rhombus with 
measured probabilities P m i,i = 1,3 (see Fig. |2}. Each 
of these distortions involves four 7—5 disclinations i.e. 
two possible dislocation- antidislocation pairs which, we 
assume, occur independently. For a free (Vo = 0) two 
dimensional system dislocation core energy E l c can be 
found through the relation 3 ^ 



o o 



II = exp(-/3 2E t c )Z{K) 



(3) 



where II = 53i=i ^mi an d Z(K) is the "internal partition 
function" incorporating all three types of degenerate ori- 
entations of dislocations, 
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where Jo is a modified Bessel function, K — 0KOq is a di- 
mensionless Young's modulus renormalized over phonon 
modes, ao being the lattice parameter and r m i n is the sep- 
aration between dislocation-antidislocation above which 
one counts the pairs. The above expression for Z(K) 
and Eq.@ have been used previously in simulationsiiSi 
of phase transitions of 2-d systems in absence of any ex- 
ternal potential to find the dislocation core energy E\. 

The probabilities for occurrence of the dislocation pairs 
of a specific type themselves P^i (Fig. |3J| which are pro- 
portional to the square of the fugacities, can be com- 
puted easily. The probability of dislocation pairs of type 
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FIG. 2: This diagram depicts all the possible dislocation 
generating moves that we reject. Starting from the triangular 
lattice shown in the centre (the dotted lines show the potential 
minima), in all, there can be three types of dislocation- pair 
generating moves shown as A, B & C. The numbers 7 and 5 
denote the positions of two types of disclinations having seven 
nearest neighbours and five nearest neighbours respectively. 
Only those bonds, which are necessary to show distortions due 
to the generation of disclination quartets, have been drawn. 
The rhombi near each of the distorted lattice denote the unit 
cells and open arrows from 7 — > 5 show the direction of dislo- 
cation generating moves. The probabilities of these moves are 
fmi(A), P m 2(B) and P m 3(C). Corresponding Burger's vectors 
(filled arrows) are bisectors pointing towards a direction ro- 
tated counter-clockwise starting from 7 — >■ 5 directions and 
are parallel to one of the lattice planes. Notice, the sepa- 
ration between Burger's vectors of a pair along the glide di- 
rection (parallel to the Burger's vectors) is a single lattice 
separation (oo) and within this construction it is impossible 
to draw Burger's loop that can generate non-zero Burger's 
vector. Depending on which of the two possible disclination 
pairs separate out any one dislocation- antidislocation pair 
will be formed. 



I is P dl = ±(P m2 + P m3 - P m i) and that of type II is 
Pd2 = \{P m i + P m 2 + Pm3 - 2P A1 ) = P ml /2. The va- 
lidity of the above expressions can be clearly seen from 
FigE 

An argument following the lines of Fisher et. alM shows 
that the dislocation probability (number density of dis- 
location pair per unit cell) for our system, 



P d i = exp(-/3 2E c )Z(K xy ) 



(4) 



where 2E C is the core energy and Z{K xy ) is the inter- 
nal partition function of dislocation pair of type I (single 
orientation) . 
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with K xy — [3K xy and A c — s/Zoq/2 being the area of 
an unit cell in the undistorted lattice. We choose r m i n = 
2dQ. At this point this choice is arbitrary. We give the 
detailed reasoning for this choice at the end of section 
IIIII EqQ] and Eq[5] straightaway yield the required core 
energy E c . The corresponding fugacity contribution to 
RG flow equations (Eq|2J is given via 

y' = An^P dl /Z{k xy ) (6) 

In the above, the following assumption is, however, im- 
plicit. Once a nearest neighbor bond breaks and a poten- 
tial dislocation pair is formed, they separate with proba- 
bility one^. This assumption goes into the identity Eq0] 
as well as in Eq[2f . Taking the rejection ratios due to 
bond- breaking as the dislocation probabilities, as well, 
require this assumption^. 

The same restricted Monte Carlo simulation can be used 
to find out the stress tensor, and the elastic moduli from 
the stress-strain curves. The dimensionless stress tensor 
for a free (Vq = 0) system is given by2£, 

(- ^{^ r ¥) + < ? » 

where i, j are particle indices and A, v denote directions 
x, y; (/>(r lJ ) is the two- body interaction, S/d 2 is the 
dimensionless area of the simulation bojsS. 



III. RESULTS AND DISCUSSION 

In this section we present the calculation of the phase 
diagram for three different 2-d systems, namely hard 
disks, soft disks and a system of colloidal particles 
interacting via the DLVO (Derjaguin-Landau-Verwey- 
Overbeek) 3 -' 4 ? potentials. We discuss, first, the calcula- 
tion of the phase diagram for a two dimensional system 
of hard disks, in detail. The bulk system of hard disks 
where particles i and j, in 2-d, interact via the potential 
0(r y ) = for r*- 7 > d and </>(r y ) = oo for r y < d, where 
d is the hard disk diameter and r y = |r J ' — r l | the relative 
separation of the particles, is known to melt^SdLiSi^ 
from a high density triangular lattice to an isotropic liq- 
uid with a narrow intervening hexatic phaseiiLiSiS. The 
hard disk free energy is entirely cntropic in origin and the 
only thermodynamically relevant variable is the number 
density p = N/V or the packing fraction rj = (n/A)pA 2 . 
Simulations 32 , experimental^ and theoretical studies 
of hard disks show that for r\ > .715 the system exists 
as a triangular lattice which transforms to a liquid be- 
low r] = .706. The small intervening region contains a 
hexatic phase predicted by the KTHNY theorjiiiiSi of 
2-d melting. Apart from being easily accessible to the- 
oretical treatments^, experimental systems with nearly 
"hard" interactions viz. sterically stabilized colloids^ 
are available. 
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FIG. 3: In this plot the O symbols correspond to Pdi, the 
probability for type I dislocations and the X symbols to Pd2 
the probability for type II dislocations obtained from the P m i 
(see text and Fig|5J for various r\ values, arrow denoting the 
direction of increasing r\{= .69, .696, .7029, .71). Pdi for i = 
1,2 are expressed in units of 10 4 . These probabilities are 
plotted against the potential strength /3Vq. Note that for 
PVo > 1, the probability for type I dislocations is larger than 
that of type II. The dots and solid lines are only guides to 
eye. 



In presence of a periodic external potential, the only 
other energy scale present in the system is the relative 
potential^ strength (3Vq. If the modulating potential is 
commensurate with the spacing between close- packed 
lines, the elastic free energy of this system in it's solid 
phase follows Eq^ and the corresponding renormaliza- 
tion flow equations are given by Eq|21 

We obtain the bare y' and x' from Monte Carlo simula- 
tions of 43 x 50 = 2150 hard disks and use them as initial 
values for the numerical solution of Eqs. J3J). The Monte 
Carlo simulations for hard disks is done in the usual^i 
way viz. we perform individual random moves of hard 
disks after checking for overlaps with neighbours. When 
a move is about to be accepted, however, we look for the 
possibility of bond breaking as described in the previous 
section (Fig|2J). We reject any such move and the rejec- 
tion ratios for specific types of bond breaking moves give 
us the dislocation probabilities of type I and II, sepa- 
rately (Fig|2Jl. From Fig[3]it is clear that the probability 
of type II dislocations i. e. P$i drops down to zero for all 
packing fractions at higher potential strengths (3Vq. The 
external potential suppresses formation of this kind of 
dislocations. For small /3Vo on the other hand, the prob- 
abilities of type I and type II dislocations are roughly the 
same. This should be a cause of concern since we neglect 
the contribution of type II dislocations for all 0Vq. We 
comment on this issue later in this section. 

Using Eq|H| and EqJSl along with the identity r m i n = 2ao 
gives us the initial value y' to be used in renormalization 
flow Eq[21 provided we know K xy . Again K xy gives x' 
straightaway. To obtain that we need to calculate the 
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FIG. 4: Plot of f3d 2 a xx vs. 8/d at a strain value e xx = .02 for 
packing fraction r\ — .7029 and potential strength Vo = 1. A 
second order polynomial fit (solid line) gives lim^^o f3d 2 a xx = 
-6.23 . 
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FIG. 5: A typical stress-strain curve used to obtain the Young 
modulus from a linear fit (solid line). The graph is plotted 
at rj = .7029, Vq = 1.0. The fitted Young's modulus f3Kd 2 = 
54.84. 



Young modulus K and shear modulus /i. 

Let us again go back to Eq0 the expression for stress 
tensor. For hard disk potentials the derivative d(f)/dr 1 ^ 
becomes a Dirac delta function and the expression for 
stress can be recast into^l 

( 8 ) 

The presence of Dirac delta function 8(r 1 ^ — d) in the 
above expression requires that the terms under the sum- 
mation contribute, strictly, when two hard disks touch 
each other i.e. r lJ = r = a. In practice, we imple- 
ment this, by adding the terms under summation when 
each pair of hard disks come within a small separation 
r = a + 8. We then find (3a\ u d 2 as function of 8 and fit 
the curve to a second order polynomial. Extrapolating to 
the 8 — > limit obtains the value of a given component 
of stress tensor at each strain value e\^L. 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



5/d 

FIG. 6: Plot of (3d 2 a xy vs. 8/d at strain value e xy = .08 at the 
packing fraction r\ — .7029 and potential strength Vo = 1. A 
second order polynomial fit (solid line) gives lima^o Pd 2 a xy = 
1.092 . 
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FIG. 7: A typical stress-strain curve used to obtain shear 
modulus from a linear fit (solid line) . The graph is plotted at 
r) = .7029, Vb = 1.0. The fitted shear modulus fifid 2 = 13.53. 



For completeness, now we show how we calculate the 
two relevant stress-tensors : a xx at a given longitudinal 
strain e xx in Fig. 0] and a xy for a shear strain e xy in 
Fig.0 We thus calculate the stress at each value of strain 
and from the slopes of stress-strain curves find out the 
bare Young-modulus (3Kd 2 (Fig. [^J and shear-modulus 
f3 fid 2 (Fig. 0). To obtain the relevant elastic moduli we 
give first an elongational strain in x- direction which is 
parallel to the direction of potential minima to obtain K 
and subsequently a shear in the same direction to obtain 
(i. Any strain that forces the system to ride potential 
hills will give rise to massive displacement modes which 
do not contribute to elastic theory. 

From these elastic moduli we get the 'bare' K xy (and 
hence x' — ttK xv — 2, see section [nj. This is also re- 
quired to complete the computation of y' . In Fig. |S| 
we have plotted x' and y' the bare values of x' and y' 
for various potential strengths f3Vo at packing fraction 
r/ = .7029 along with the separatrices for the linearized 
and the non- linear flow equations (Eq. [21 . The line of 
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FIG. 8: The initial values of x' and y' obtained from the 
elastic moduli and dislocation probability at n — .7029 plot- 
ted in x' — y' plane. The line connecting the points is a 
guide to eye. The arrow shows the direction of increase in 
PV {= .01, .04, .1, .4, 1, 4, 10, 40, 100). The dotted line denotes 
the separatrix (y' = x') incorporating only the leading order 
term in KT flow equations. The solid curve is the separa- 
trix when next to leading order terms are included. In I — > oo 
limit any initial value below the separatrix flows to y' = line 
whereas that above the separatrix flows to y — *• oo. The inter- 
section points of the line of initial values with the separatrix 
gives the phase transition points. The plot shows a freezing 
transition at /3Vo = .1 followed by a melting at /3Vo = 30. 



initial conditions is seen to cross the non-linear separa- 
trix twice (signifying re-entrant behaviour) while cross- 
ing the corresponding linearized separatrix only once at 
high potential strengths. The phase diagram (Fig. I10|) 
is obtained by computing the points at which the line 
of initial conditions cut the non-linear separatrix using 
a simple interpolation scheme. It is interesting to note 
that within a linear theory the KT flow equations fail 
to predict a RLIF transition. Performing the same cal- 
culation for different packing fractions rj we find out the 
whole phase diagram of RLIF in the rj- (3Vo plane. 

The numerical errors in the phase diagram are calculated 
as follows. The quantity (3K xy varies linearly with rj at 
all potential strengths. Therefore the numerical error in 
r) is proportional to the error in /3K xy (see Fig. |5J). Using 
all these we obtain the RLIF phase diagram for hard disk 
systems (Fig lTUfr . 

Further, comparing with previous computations^^ of 
the phase diagram for this system (also shown in Fig. llO|) 
we find that, within error- bars, our results agree at all 
values of r\ and (3Vo with the results of W. Strepp et. 
al~. Whereas, in numerical details, they disagree with 
the results of C. Das et. alm^, though even these re- 
sults show RLIF and are in qualitative agreement with 
ours. This validates both our method and the quanti- 
tative predictions of Ref»2&. The effect of higher order 
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FIG. 9: For hard disk system f3K xy varies linearly with rj. 
Data plotted at Vo = 1. The solid line is a linear fit to the 
form f(x) = a + bx with a = -7.37 and b = 11.76. At 
each Vo the error in K X y determines the error in n: Sn/n = 
\l + a/r,b\{5K xy /K xy ). 
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FIG. 10: The phase diagram of the hard disk system in the 
presence of a 1-d, commensurate, periodic potential in the 
packing fraction (rj) - potential strength (/3Vo) plane. The 
points denoted by □ correspond to our RG calculation using 
the techniques described in this paper. The points denoted by 
OS and *^ are taken from earlier simulations. The vertical 
bars denote estimate of error. Our data clearly matches with 
Ref[7].The horizontal line at n — .706 denotes the calculated 
asymptotic phase transition point at (3Vo = oo. 



terms in determining non-universal quantities has been 
pointed out earlier— for the planar rotor model but in the 
present case their inclusion appears to be crucial. Nev- 
ertheless, we expect our procedure to break down in the 
(3Vo —> limit where effects due to the cross-over from 
a KT to a KTHNYii-iS transition at f3V = become 
significant. Indeed, as is evident from Fig.|31for f3Vo < 1 
the dislocation probabilities of both type I and type II 
dislocations are similar— and the assumptions of FNR 
theory and our process (which involves only type I dis- 
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FIG. 11: Phase diagram for soft disks: □ denote our calcu- 
lation, O indicate earlier simulation data-£<—i. The vertical 
lines are the error- bars. 
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FIG. 12: Phase diagram for particles interacting via the 
DLVO potential. □ denote our calculation, O show the ear- 
lier simulation data™. The vertical lines are the error- bars. 
Error bars in our calculation being smaller than the symbol 
size are not shown. 



locations) need not be valid at small potential strengths. 
This fact is also supported by Ref^ where it was shown 
that though at /3Vq = 1000 the scaling of susceptibility 
and order parameter cumulants gave good data collapse 
with values of critical exponents close to FNR predic- 
tions, at PVq — .5, on the other hand, ordinary critical 
scaling gave better data collapse than the KT scaling 
form, perhaps due to the above mentioned crossover ef- 
fects. In the asymptotic limit of (3V$ — * oo the system 
freezes above r\ = .706 which was determined from a sep- 
arate simulation in that limit. This number is very close 
to the earlier value rj ~ .71 quoted in RefiS. As expected, 
the freezing density in the (3Vq — > 00 limit is lower than 
the value without the periodic potential i.e. r\ ~ .715. 

In a similar fashion it is possible to find out phase di- 
agrams of any 2-d system in presence of external mod- 
ulating potential commensurate with the density of the 



system. We illustrate this by calculating similar phase 
diagrams for two other systems, viz. soft disks and the 
DLVO system. The soft disks interact via the potential : 



0(r) 



1 



where r denotes the separation between particles. In 
simulations, the cutoff distance is chosen to be r c = 2 
above which the particles are assumed to be non- in- 
teracting. Apart from the external potential strength 
(3Vq the relevant thermodynamic quantity is the number 
density p — N/L x L y . In finding 'bare' elastic moduli 
from restricted simulations the stress is calculated from 
the Eq[7| As this expression does not involve any Dirac 
delta functions unlike the case of hard disks, we do not 
require any fitting like Figs 14161 to obtain the stresses. 
The elastic moduli are again found from stress- strain 
curves like Figs l5l7l The dislocation fugacity of type I is 
calculated from rejection ratio of dislocation generating 
moves. All these, at a given p value generate the initial 
conditions x' and y' in RG flow diagrams. The crossing 
of these initial conditions with the separatrix found from 
Eq0 gives the phase transition points. The phase dia- 
gram is plotted and compared with phase diagram from 
earlier simulations^^ in Figllll The error bar in p is 
found from the error in K xy , as K xy varies linearly with 
p, through the relation Sp/p = |1 + a/ pb\(SK xy / K xy ). 
The quantities a and b are found from linear fitting (of 
form a + bx) of K xy vs. p curve, at any given /3Vq. 
The phase diagram (Fig. 111( 1 again clearly shows re- 
entrance (RLIF). This is in qualitative agreement with 
earlier simulations^^ (see Fig llll . 

For charge stabilized colloids the inter-particle poten- 
tial that operates is approximately given by the DLVO 
potentials^: 

. . (Z*e) 2 ( exp(.5fcd)\ 2 exp{—nr) 



47reoe r \ 1 + .5nd 

where r is the separation between two particles, d is the 
diameter of the colloids, k is the inverse Debye screen- 
ing length, Z* is the amount of effective surface charge 
and e r is the dielectric constant of the water in which the 
colloids are floating. In order to remain close to exper- 
imental situations and to be able to compare our phase 
diagram with the simulations of Strepp et. al^ we use 
T = 293.15^, d= 1.07/jm, Z* = 7800, e r = 78. In 
experiments, the dimcnsionless inverse Debye screening 
length na s can be varied either by changing k through 
the change in counter-ion concentration or by changing a s 
by varying density^. We perform the restricted Monte- 
Carlo simulation as described in section ^ In simula- 
tions we vary k and keep the particle spacing in ideal 
lattice a s = 2.52578/j,to (density) fixed. Further, we use 
a cut- off radius r c such that, <fi(r > r c ) — where r c 
is found from the condition j34>(r c ) = .001. We find 
out phase transition points (in na s ) at different exter- 
nal potential strengths f3Vo in the same fashion as de- 
scribed earlier. The bare renormalizable quantities x' 
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FIG. 13: The figures a - d which we have drawn using the 
applet "voroglide" 50 show four steps of separation of a type I 
dislocation pair, from a separation of ao to 4ao. The shaded 
regions show the 5 — 7 disclination pairs constituting the dis- 
locations. Burger's circuits are shown in a - c. Note that 
for separations > 2ao separate Burger's circuits around each 
disclination pair give rise to non-zero Burger's vectors, giv- 
ing the dislocations their individual identity. This shows that 
the minimum meaningful separation between dislocation cores 

^mi)i — 2ttQ . 
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FIG. 14: Similar to Fig|H| The initial conditions x' Q and y' 
are plotted as a function of /3Vo- The different data sets are 
created for different values of r min . The symbols mean the 
following : O denotes data for r min = ao, □ denotes that for 
rmin = 2ao and A denote data for r min — 3ao. The dotted 
line denotes the non- linear separatrix. 



and y' are found from restricted Monte- Carlo simula- 
tions for various (3Vq at each na s . The phase transi- 
tion points are calculated from the intersection of these 
initial conditions with the separatrix found from Eq|5J 
Thus we obtain the phase diagram in Ka s — (3Vq plane 
(Fig. I12fl . (3K xy varies linearly with na s and the error 
in K xy generates the error in na s fFig ll2li through the 
relation £(rea s )/(/ca s ) = |1 + a/bKa s \(5K xy / K xy ). The 
quantities a and b are found from linear fitting (of form 
a + bx) of K xy vs. na s curve, at any given f3Vo. Though 
there is a quantitative mismatch between our data and 
that of Strepp et. alm^, our data shows a clear region in 
na s (between 15.1 and 15.2) where we obtain re-entrance 
(RLIF). This is in contrast to the simulated phase dia- 
gram of C. Das et. alm^, where they observe absence of 
re-entrance at high field strengths. We do not plot their 
data as the parameters these authors used are not the 
same as the ones used in Figll2l 

It is interesting to note that, with increase in range of 
two- body interaction potentials the depths of re-entrnace 
(in rj 1 p or na s ) decreases. This is again in agreement with 
the understanding that, the re-entrant melting comes 
about due to decoupling of the 1-d trapped layers of par- 
ticles that reduces the effective dimensionality thereby 
increasing fluctuations. With the increase in range of the 
interacting potentials this decoupling gets more and more 
suppressed, thereby reducing the region of re-entrance. 

One aspect of our study which stands out is the excep- 
tionally better agreement of our results with previous 
simulations for hard disks as opposed to systems with 
soft potentials like the soft disks and the DLVO. This 



could, in principle, be due either (a) to the failure of the 
RG equations used by us or some other assumptions in 
our calculations (b) or to unaccounted finite size effects 
in earlier simulations. While it is difficult to estimate 
the effect of (a) since RG equations to higher orders in y 
are unknown at present, we may be able to motivate an 
estimation for (b). In order to explain the discrepancy in 
the positions of the phase boundaries, we need to go into 
some details of how the phase diagrams were obtained in 
the earlier simulations. In these simulations^Si 2 ^ 6 . the 
phase boundaries were obtained from the crossing of the 
order parameter cumulants^S^i for various coarse grain- 
ing sizes. The system sizes simulated in these studies are 
the same (N = 1024). However, the range of interac- 
tion differs. To obtain an objective measure we define 
the range of the potentials £ as that at which the inter- 
action potential <j> is only 1% of its value at the lattice 
parameter. In units of lattice parameter, we obtain, for 
soft disks £ = 1.47 and for the DLVO potential £ = 1.29 
at typical screening of na s = 15. By definition, for hard 
disks £ = 1. The particles within the range of the po- 
tential are highly correlated and we calculate the num- 
ber N corr of such independent bare uncorrelated parti- 
cles within the full system size. N corr takes the values 
N corr = 1024, 473.88, 615.35 for hard disks, soft disks 
and the DLVO potential respectively. Since the effective 
system sizes are smaller for the soft potentials, finite size 
effects are expected to be larger. In this connection, it is 
of interest to note that in the same publications 2 ^ 4 * 2 ^ 6 , 
a systematic finite size analysis showed that the phase di- 
agrams shift towards higher (lower) density (kappa) for 
hard and soft disks (DLVO). A look at FigfTTl and IT2I 
should convince the reader that such a shift would ac- 
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tually make the agreement with our results better. We 
emphasize here that our present restricted simulations 
are virtually free of finite size effects since the system 
does not undergo any phase transition. 

Before we end this section, we discuss the reasons behind 
the particular choice of r m i n that we made throughout 
this manuscript. In practice, it is possible to give indi- 
vidual identification to a dislocation only when Burger's 
vector separation within a pair is > 2ao (Fig ll3[) i.e. 
Train = 2ao. For r > 2ao Burger's loops can be drawn 
around each 5 — 7 disclination pair (Fig ll3l) giving rise 
to a non-zero Burger's vector. In unrestricted simula- 
tions and in experimental situations after a disclination 
quartet is formed, they get separated out and the easy 
direction of separation is the glide direction which is par- 
allel to the Burger's vector. In Fig^|we show four steps 
of separation of such a dislocation pair of type I. After 
motivating r TO j„ = 2ao we show, in Fig ll4l the three sets 
of initial values corresponding to r m j n = ao, 2ao, 3ao 
along with the non- linear separatrix at r\ = .7029 of hard 
disk system. r m i n = ao predicts the system to be in 
solid phase for any arbitrarily small amount of external 
potential and to melt at larger (3Vq . This behaviour con- 
tradicts physical expectation that the melting density at 
f3Vo = has to be larger than that at /3Vb = oo. On 
the other hand, while r m j ra = 3ao does not produce any 
unphysical prediction, it shrinks the region of re-entrancc 
in the (3Vq direction. Therefore r m i„ = 2ao is the min- 
imum value for r m i n that could be chosen to produce 
physically meaningful results and this choice remains in 
closest agreement with simulation data. 

IV. CONCLUSION 

We have presented a complete numerical renormal- 
ization group scheme to calculate phase diagrams for 



2-d systems under a commensurate modulating poten- 
tial. We have used FNR theory along with this scheme 
to calculate phase diagrams for three different systems, 
namely, the hard disks, the DLVO and the soft disks. In 
all the cases we have found laser induced freezing followed 
by a re-entrant laser induced melting. We show that the 
re-entrance behavior is built into the 'bare' quantities 
themselves. We find extremely good agreement with ear- 
lier simulation results. In particular the phase diagram 
for hard disk comes out to be exactly the same as found 
from one set of earlier simulations 23 . To obtain the cor- 
rect phase diagram, however, flow equations need to be 
correct at least upto next to leading order terms in the 
dislocation fugacity. Our results, especially for small po- 
tential strengths, is particularly sensitive to these terms. 
Cross-over effects from zero potential KTHNY melting 
transition are also substantial at small values of the po- 
tential. 
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